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Abstract 



Over the past decade, the numerical modehng of the magnetic field evolution in astrophysical scenarios has become an 
' ^ increasingly important field. In the crystallized crust of neutron stars the evolution of the magnetic field is governed by 
the Hall induction equation. In this equation the relative contribution of the two terms (Hall term and Ohmic dissipation) 
varies depending on the local conditions of temperature and magnetic field strength. This results in the transition from 
^^the purely parabolic character of the equations to the hyperbolic regime as the magnetic Reynolds number increases, 
which presents severe numerical problems. Up to now, most attempts to study this problem were based on spectral 
<H methods, but they failed in representing the transition to large magnetic Reynolds numbers. We present a new code 
based on upwind finite differences techniques that can handle situations with arbitrary low magnetic diffusivity and it is 
Suitable for studying the formation of sharp current sheets during the evolution. The code is thoroughly tested in different 
limits and used to illustrate the evolution of the crustal magnetic field in a neutron star in some representative cases. 
Our code, coupled to cooling codes, can be used to perform long-term simulations of the magneto-thermal evolution of 
neutron stars. 
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1. Introduction 

The magnetic field (MF) is a key agent to explain many 
'of the observed properties of neutron stars (NSs). In or- 
der to explain NS rotational energy losses by dipole radi- 
'ation, a simple back-of-the-envelope estimate gives that 
JSfSs must have an external dipolar component ranging 
from 10"'^'' to 10^^ G, but little is known about its inter- 
nal geometry and evolution. In the standard NS evolution 
"model, a few hours, or at most a few days after birth, 
a solid crust about 1 km thick (compared to the typical 
10 km radius) is formed. Under these conditions, ions in 
,the lattice have very restricted mobility and conduction is 
■governed by the electrons. At the large densities of a neu- 
,tron star crust (10^° — 10^^ g cm~^), matter is completely 
■pressure-ionized and electrons can flow trough the lattice. 
The equation describing the evolution of such system is the 
■generalized Hall induction equation for a one-component 
plasma (electrons). This limit is also known as Electron 
MHD (EMHD), and describes a one-component charged 
fluid which density does not vary with time. 

The evolution of the crustal MF is thus regulated by 
Ohmic dissipation, with timescales of 10^ yr, and the Hall 
term. The latter causes a transfer of energy to smaller 
scales and the coupling between toroidal and poloidal field 
components, on timescales of 10^ — 10^ yr. The Hall- 
driven evolution is probably at the origin of the observed 



strong activity in young magnetars, i.e., NSs with ultra- 
strong magnetic fields. Therefore, the interest in modeling 
the internal evolution of t he MF is continuously growing 
dPons and Geppertl. [2OO7I: IHovos et all l2008t IPons et al 



2009t IShabaltas and Lari2012t ). The MF in the crust of 



NSs also directly affects the microphysical processes that 
govern the ther mal evolution of the crust. For this reason . 



2007: Pons et al 



previous works ( Aguilera et all . 120081 : iPons and Geppert 



20091) were directed toward a fully cou 
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pled m agneto-thermal evolution of the NS crust. lAguilera et al 
( 20081 ) developed a 2D (axial symmetry) cooling code tak- 
ing into account the interplay between temperature and 
MF. However the decay of the latter was simulated by a 
phenomenological, homologous analytical formula, with- 
out solving the induction equation because of numerical 
limitations in the treatment of the Hall term. One of the 
important parameters is the resistivity, which strongly de- 
pends on the local temperature. For typical temperatures 
(a few 10^ K) in a middle-age NS of 10^ - 10^ yr, the Hafl 
term dominates over Ohmic dissipation when the MF is 
> 10^"* G. However, as the star cools down and the Ohmic 
timescale becomes very long (billions of years), the Hall 
term may dominate even for much lower values of the MF 
strength, which in turn requires a better numerical ap- 
proach. 

The numerical challenge is to be able to follow the 
evolution of such system, in which important parameters 
(density, resistivity) vary several orders of magnitude across 
the spatial domain, but also during the temporal evolution. 
A multi-purpose code must be able to work in both the 
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purely diffusive regime and in the limit of large magnetic 
Reynolds number when the Hall term dominates. It must 
also be stable to follow the evolution for many diffusion 
timescales (up to hundreds of Hall timescales). Histori- 
cally, spectral methods had been used to solve the Hall 
induction equation in si i nplifie d constant density layers 



( Hollerbach and Riidieeii 120021 ) . but such attempts were 



always restricted to a magnetic Reynolds number not ex- 
ceeding 200, since fully spectral codes systematically have 
unsurmountable problems dealing with structures where 
discontinuities or very large gradients of the variables ap- 
pear, which is a natura l consequence of the equations. 
Pons and Gepper^ ( 2007 ) presented a code solving the Hall 
induction equation in a realistic crust using an alternative 
approach (spectral in angles but finite differences in the 
radial direction). This was a significant improvement to 
previous work, but still could not work in the limit of van- 
ishing electrical resistivity. For this reason, the only long- 
term fully coupled 2D magn eto-thermal evolution simu- 
lations available up to now (IPons et al.l . l2009l) were re- 



stricted to the purely diffusive case. 

In this paper we present a new code based on upwind, 
finite difference schemes that can handle the Hall term in 
the induction equation for vanishing physical resistivity. 
Preliminary results using the code described in this paper, 
and confirming the occurrence of the Hall instability i n 
neutron stars were presented in IPons and GeppertI ()2010l ). 
The shock-capturing character of the method employed 
ensures that the formation of current sheets is properly 
modeled. The code can follow the evolution of complex 
geometries with large gradients and discontinuities over- 
coming intrinsic problems of previous studies. The paper 
is organized as follows. In Sj2]we describe the generalized 
induction equation, and discuss the importance of the re- 
sistive and Hall terms in the system. In fJ2] we present 
the numerical code, describing the grid and the scheme 
used to solve the induction equation. In 21 '^^ present the 
battery of tests we have performed, and show some rep- 
resentative models. In fj5] we consider the evolution of a 
purely toroidal field in a realistic NS crust, with particu- 
lar attention to the formation of discontinuities. In fJS]we 
present a general case with a realistic mixture of poloidal 
and toroidal MF. In SJ7] we summarize our findings and 
discuss future applications of physical interest. 

2. The EMHD limit and the Hall induction equa- 
tion 

In the EMHD approximation, there is the implicit as- 
sumption that the timescale of variation of electromagnetic 
field is much larger than the timescale of collisions inside 
plasma. Therefore we can neglect displacement currents 
in Ampere's law and the current density is simply 



J = —V X B 
An 



(1) 



In a one-component plasma (in this case, electrons), the 
velocity is not an independent variable, but instead it is 
related to the current and charge densities by 



J 



A-Ken, 



-V X B , 



(2) 



where rie is the electron number density, e the elementary 
charge, and c the speed of light. 

For a conducting fiuid, moving with velocity and in 
the presence of a MF, the electric field [E) is related to 
the current density through the generalized Ohm's law 



J = a\ E+-xB 



(3) 



where a is the electrical conductivity. 

The evolution of this syste m is only governed by the 
HaU induction equation (e.g. IPons and Geppert ( 2007 ) 
and references therein). 



{V xB)x B} (4) 



dB ^ r ^ ^ c 

which can be derived from the induction equation 

- = -cVxi?, (5) 
where, from Eqs. ([2]) and (O, the electric field is given by 



cE ^7]^ X B 



Airen, 



■(V X B)x B. 



(6) 



Since the magnetic diffusivity, 77 = /Aira, depends on the 
temperature, the magnetic and thermal evolution in the 
crust of a NS are coupled. The first term on the right-hand 
side of Eq. (g]) is responsible for the Ohmic decay. The 
second term is the Hall term, and its non-linear character 
is in the root of all difficulties one finds when numerical 
methods are used to solve the Hall induction equation. 

To compare the relative importance of the two terms 
(Ohmic dissipation and Hall), we can write the Hall in- 
duction equation @ as follows: 



dB 
'dt 



= -V x 



(ryjv X B + 7^„(V X B) X (7) 



where b is the unit vector in the direction of the MF and 
TZm = cB / Airerieri is the magnetic Reynolds number. TZm 
is an indicator of the relative importance between the two 
terms on the right-hand side. 

We define the characteristic Ohmic dissipation timescale 

LVv , (8) 



and Hall timescale 



Th = AneUeL'^/cB = Td/TZr. 



(9) 



where L is a typical length scale. With these definitions, 
TZm is the ratio between the dissipation and Hall timescales. 
It also corresponds to the so-called "magnetization param- 
eter" (wBTe), usually found in the literature of plasma 
physics, where ujb is the cyclotron frequency and Te the 
electron relaxation time. 



2 



3. The numerical code 

In order to understand the dynamical evolution of the 
system and to design a successful numerical algorithm, it 
is important to know the mathematical character of the 
equations and to identify the wave modes. The magnetic 
Reynolds number determines the relative importance of 
the Hall term and Ohmic dissipation, and defines the tran- 
sition from a purely parabolic equation {TZm <C 1), to a hy- 
perbolic regime {TZm ^1)- In the case of Hall-dominated 
evolution, the problem becomes a "doubly" constrained 
multidimensional advection equation because, in addition 
to the W-B — constraint, the velocity field is proportional 
to the current (Eq. [5]), i.e., to the derivatives of the field. 
This implies non-linearity and coupling between poloidal 
and toroidal components. An initially pure toroidal field 
will not develop a poloidal component, but any poloidal 
field will induce the formation of toroidal components. 
When both components are present, magnetic energy can 
be transferred between them. Although the Hall term it- 
self conserves energy, the plausible creation of small-scale 
structures may accelerate the Ohmic dissipation. 

The Hall term also in troduces two wave modes into 
the system. Huba (j2003h has shown that, in a constant 
density medium, the only modes of the Hall EMHD equa- 
tion are the whistler or helicon waves, which are transverse 
field perturbations propagating along the field lines. In 
presence of a charge density gradient, additional Hall drift 
waves appear. These are transverse modes that propa- 
gate in the B x Vrie direction. It has also been noted 
that the presence of charge density gradie nts results in a 



Burgers-like term ( Vainshtein et al. ■ I2OOOI) . Furthermore 



even in the constant density case but without planar sym- 
metry, the evolution of the toroidal component of the MF 
also contains a quadratic term that re sembles the Burg- 
ers equation (jPons and Gepoertl . 120071 ) with a coefficient 
dependent on the distance to the axis. A natural out- 
come of the presence of Burgers-like terms is the forma- 
tion of "shocks" (current sheets). The evolution for the 
poloidal field is more complicated, as it includes higher 
order derivatives in the non-linear terms. All the above 
issues must be taken into account when an efficient and 
robust numerical method is designed to solve the Hall in- 
duction equation in the limit of strong magnetization. To 
overcome all previous limitations found when techniques 
based on spectral methods have been used, we decided 
to use well-known upwind, conservative schemes. These 
methods are of very general use, from problems involving 
the simple ID Burgers equation to high resolution shock 
captu ring schemes very successfully used in MHD prob- 
lems ([Anton et aT . 2006 : Giacomazzo and RezzollaL 2007 : 



Cerda-Duran et al.l . l2008h . 



3.1. The numerical grid: staggered grids 

Staggered grids (|Yeel . ll96^ are c ommonly used with fi- 
nite difference time domain methods ([Taflove and Brodwin , 




Figure 1: Location of the variables on a staggered grid in spherical 
coordinates and for the axisymmetric case. Solid lines delimit the 
edges of the surface S^,. 



1975il to solve the Maxwell's equations. The MF compo- 
nents are defined at the center, and normal to each face of 
the cell, while the electric field is defined along the edges 
of the cell. With these definitions, the W ■ B — condi- 
tion can be satisfied to machine error (see below). The 
components of V x B, as the electric field, are also defined 
along the edges of the cell. Stokes' theorem allows us to 
write the curl components in terms of the MF components 
defined at the cell faces. 

As an example, we show in Fig. [T] the location of the 
variables in a numerical cell in spherical coordinates (r, 9, (p) 
and assuming axial symmetry. In this case our grid maps a 
meridional section of the star. In this paper we will restrict 
to 2D simulations, but we will present tests and results in 
both Cartesian and spherical coordinates. 

3.2. Conservation form and divergence-preserving schemes. 

One important feature of numerical schemes designed 
for hyperbolic problems is the formulation of the equations 
in conservative form. This ensures a correct propagation of 
waves and the preservation (by construction) of the diver- 
gence constraint. Integrating over the surface of a cell face, 
normal to the a direction, and applying Stokes theorem to 
Eq. ([S]), the space-discretized equation for the magnetic 
flux <I>Q, can be written as 

1 



dt 



fe 



Ed 



(10) 



The circulation of the E-iie\d around a face boundary has 
been approximated by the sum J^k^kh, where Ek are 
the central values of the electric field on the edges, Ik is 
the length of the edge, and k runs over the four edges of 
the face. The normal components of the MF to a given 
face are assumed to be average values over the surface, 
Ba — ^af^a, wherc is the face area. 

Making use of Gauss' theorem, the numerical diver- 
gence can be evaluated, for each cell with volume AV, as 
follows: 

1 



V • B 



AV ^ 



(11) 



where the plus and minus superscripts indicate the fluxes 
through the two opposite faces in a given direction a (in 
the 2D example of Fig. [1] a = r,6). With this defini- 
tion, the divergence-preserving character of the methods 
using the conservation form and advancing in time mag- 
netic fluxes, instead of MF components, becomes evident: 
taking the time derivative of Eq. (fTT|) . and using Eq. (fTO)) . 
every edge contributes twice with a different sign and can- 
cels out. 

3.3. Cell reconstruction 

The original Godunov's method is well known for its 
ability to capture discontinuous solutions, but it is only 
first-order accurate. This method can be easily extended 
to give second-order spatial accuracy on smooth solutions, 
but still avoiding non-physical oscillations near disconti- 
nuities. To achieve higher order accuracy, we can use 
a reconstruction procedure that improves the piecewise 
constant approximation. The simplest choice is a piece- 
wise linear function in each cell. A very popular choice 
for the slopes of the linear reconstructed function is the 
monoton i zed ce ntral- difference limiter (MC), proposed by 
van Leeii (|l977n . Given three consecutive points Xi^i, Xi, Xi^i 
on a numerical grid, and the numerical values of the func- 
tion fi-i, fi, fi+i, the reconstructed function within the 
cell i is given by f{x) = J{xi) -I- a{x — Xi), where the slope 
is 



a = mmmoc 



fi+l — fi- 



'1 2-^^+1 



.fi 2 ■f^ — 1 



•Xi—\ XiJ^\ Xi Xi X^—X 

The minmod function of three arguments is defined by 

min(a, b, c) if a,b,c > 0; 



minmod(a, 6, c) = 



max(a, &, c) ifa, 6, c<0; (12) 
otherwise. 



We reconstruct the MF circulation elements, Ck — B^lk- 
From this, we directly obtain the components of (V x B) by 
means of Stokes' theorem and, when needed, we recover Bx 
dividing the reconstructed circulation by the local length 
element. 

From Eqs. dH) and ([S]), the electric field is: 



47r 



1 



E = ^riJ ~ -vx B (13) 

Current and electric field components are always defined 
at the same location, but the Hall term includes prod- 
ucts of tangential components of the MF and velocity not 
always defined at the same edges. For such terms, we 
evaluate the needed components of the velocity by linear 
interpolation of the closest neighbor values and we take 
the upwind components B^ of the reconstructed value of 
MF at each interface. For example, in the axisymmetric 
case and considering the evolution of the poloidal compo- 
nents, the contributions of Er and Eg to the circulation 
cancel out and we only need to evaluate the contribution 
of E^p, which is calculated by 



47r 



-{vrBg^^ -veBr) 
c 



(14) 




Figure 2: Illustration of the procedure of calculation of the electric 
field: location of the components of velocity (red arrows) and MF 
(blue) involved in the definition of the Hall term of Ey, (black dot). 



In Fig.[5]we explicitly show the location of E^ (black point) 
and the location on the staggered grid of the quantities 
needed for its evaluation. The average values of Vr and vg 
are calculated taking the average of the two closest neigh- 
bors. In the example, both Vr and vg are positive, there- 
fore the upwind values of i?^' and are the reconstructed 
values from the bottom and left sides, respectively. 

3.4. Time advance and Courant condition 

We use an explicit, first-order time advance method 
with some intermediate corrections to improve the stability 
of the scheme. In explicit algorithms, the time step is lim- 
ited by the Courant condition, that avoids that any wave 
travels more than one cell length on each time step. Since 
we want to evolve the system on long (Ohmic) timescales, 
the Courant condition is particularly restrictive when TZm S> 
1. At each time step, we estimate the Courant time (ic) 

by 



tr = min 



c|VxS|^ 



(15) 



where Al is the minimum length of the cell edges in any 
direction and {i,j) denote the numerical cells. We use a 
time step 



(16) 



where kc is a factor < 1 . The time step can vary by orders 
of magnitude during the evolution, becoming very small 
when the Hall term dominates and, in particular, where 
locally strong current sheets are formed. In a typical sit- 
uation, the Ohmic timescale is of the order of 10^ yr, the 



Hall timescale is 10^ — 10^ yr, and the Courant condition 
limits the time step between 10^'^ and 1 yr for a typical 
number of grid points of 50 x 50. 

After calculating all the line integrals of E along all 
edges, the time advance proceeds in two different steps. 
Hereafter we denote by the subscript t the component of 
the MF in the symmetry direction (the toroidal component 
in spherical coordinates, with axial symmetry), and by the 
subscript p the other two components (poloidal). First we 
advance the t component of the MF, with a particular 
treatment of the quadratic term in Bt (see S I4.1.3p . With 
the updated value of this component, the electric field com- 
ponents are recalculated and then used to advance the p 
components of the MF. Schematically, the sequence of the 
time advance from t„ to tn+i = tn + At is the following: 

• starting from all currents and electric field com- 
ponents are calculated 

^ J" ^ E"; 

• is updated: ^ 5^"+^; 

• the new values -B""*"^ are used to calculate the mod- 
ified current components and Et- 



B 



n+l 



T* 



• finally, we use the values of E^ to update the remain- 
ing MF components 



This two-step advance favors the stability of the method, 
as already p ointed out for a 3D prob l em in C artesian co- 
ordina tes by lO'SuUivan and Downed ( 20061 ). iToth et al 



further discussed that the two-stage formulation 
is equivalent to introduce a fourth-order hyper-resistivity 
term. In our case, since the t-component is advanced ex- 
plicitly, the hyper-resistive correction only acts on the evo- 
lution of the p-components. The correction introduced by 
the intermediate step is: 



<5B;^+i = -cAtV X SEt 



(17) 



where 6Et 
that 



6b;^ 



Et — After some algebra, one can show 



(cAt)2v X 



■{[Vx(Vx{ry(Vxi3r) 



iTTeUf 



[(V X S") X - (V X Bt) X B't]})] X Bn 



From this expression, we can see that the additional correc- 
tion given by 6Et contains third and fourth-order spatial 
derivatives and scales with (At)^, which is characteristic 
of hyper-resistive terms. We have found a significant im- 
provement in the stability of the method when comparing 
a fully explicit algorithm with the two-steps method. The 
method is always stable, if a sufficiently small value of fee 
is used (typically 10~^-10~^). 



3.5. Energy balance 

The magnetic energy balance equation for Hall EMHD 
can be expressed as: 



(18) 



where Qj = ATiri,P /c? is the Joule dissipation rate and S = 
cE X B/Att is the Foynting vector. During the evolution, 
the magnetic energy in a cell can only vary due to local 
Ohmic dissipation and by interchange between neighbor 
cells (Foynting flux). 

Integrating Eq. (fTS|) over the whole volume of the nu- 
merical domain, we obtain the balance between the time 
variation of the total magnetic energy £}, = Jy{B^ /87T)dV, 
the total Joule dissipation rate Qtot = Jy QjdV, and the 



Foynting flux through the boundaries St, 



JdV 



S-dY.: 



dt 



Eh 



St„t = 



(19) 



A necessary test for any numerical code is to check the 
instantaneous (local and global) energy balance given by 
Eqs. (IT^ and ([TO)) . Any type of numerical instability usu- 
ally results in the violation of the energy conservation. 
Therefore a careful monitoring of the energy balance is a 
powerful diagnostic test. 

4. Numerical tests. 

We have performed a battery of tests, with special at- 
tention to the conservation of the total energy and the 
numerical viscosity of the method employed. In order to 
understand the different aspects of Hall MHD, in this sec- 
tion we show some illustrative examples of 2D Cartesian 
tests, reproducing whistler and Hall drift waves. Finally, 
we will discuss the decay of Ohmic modes in a constant 
density sphere to compare against the analytic solutions 
for the purely resistive {TZm = 0) case. 

4.1. 2D tests in Cartesian coordinates. 
4.1.1. Whistler waves. 

The first test we have performed is to follow the cor- 
rect propagation of whistler (or helicon) waves in a two- 
dimensional slab, extending from z = — L to z = +L in 
the vertical direction. We impose periodic boundary con- 
ditions in the x-direction. All variables are independent of 
the j/-coordinate. In the case of constant density, Ug = Uq, 
and T] = (i.e. infinite magnetic Reynolds number), the 
only modes present in the system are whistler waves. We 
consider the Hall induction equation for the following ini- 
tial MF 

Bx — Bq + BiCos{kxz)cos{kxx); 

By = V2Bisiii{kxz) cos{kxx); 

Bz = Bism{kxz) sm{kxx); (20) 
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Figure 3: Evolution of the initial configuration defined by Eqs. II20II 
with Bq = 10^ Bi and kxL = tt at three different times (in units 
of tq). Arrows show the perturbed Bx (subtracting Bq), and B^- 
components, while the color scale represents the By component 
(red/black positive, yellow/white negative). 



where — mr/L, n = 1,2, and Bi <C Bq. We define 
the reference Hall timescale as 



To = 



cBn 



(21) 



This case admits a pure wave solution confined in the ver- 
tical direction and traveling in the a:-direction with speed 



-V2 k^Bo = -V2 



L^kr 



. .- , (22) 

47reno tq 

The evolution of the initial configuration defined by 
Eqs. (PO)) with Bq — 10'^ Bi, and k^L = tt at three different 
times (in units of tq) is shown in Fig. [31 and for a 200 x 50 
grid. The travel time to cross over the whole domain is 
t = 0.9 Tq, thus the perturbations have crossed through the 
horizontal domain several times, without apparent dissipa- 
tion or shape change. The code runs for hundreds of Hall 
timescales without any indication of instabilities, despite 
the fact that electrical resistivity is set to zero. The mea- 
sured speed of the whistler waves is 0.9992 the analytical 
value for this particular resolution. We have also checked 
that the correct scaling of the propagation velocity (linear 
with kx and Bq) is recovered. 



Figure 4: Horizontal section (z = 0) of the evolution of the initial 
configuration defined by Eq. I I23I I with Bq = 10^ Bi, L = 1 and 
kxL = 7r/2 at t = 0,22, and 44 (in units of tq). The perturbation 
crosses over the entire domain in 20 tq. 



4.1.2. Hall drift waves. 

As a second test in Cartesian coordinates, we set up 
the following initial configuration: 

= 0; 

By = Bq + BiCOs{kxx); 
B, = 0; 

on a stratified background in the z-direction with 



(23) 



HQ 



l + fiz 



(24) 



where ng is a reference density to which we associate the 
same Hall timescale tq defined in Eq. (|2T1) . z G [— L, L], 
and /? is a parameter with dimensions of inverse length. 
We impose again periodic boundary conditions in the x- 
direction. In the z-direction, we copy the values of the 
MF in the first and last cells {z = =fL) on the first neigh- 
bor ghost cell, to simulate an infinite domain. For small 
perturbations {Bi <^ Bq), this configuration must induce 
Hall drift waves traveling in the x-direction with speed 



Vhd 



Tq 



(25) 



For the particular model shown in Fig. 21 with Bq = 
10'^ Bi, kxL = 7r/2, and j3L — 0.2, this gives a horizontal 
drift velocity of Vhd = — 0.2_L/ro. We show the evolution 
of By at three different times, when the perturbation has 
crossed over the numerical domain 1.1 and 2.2 times, re- 
spectively. The resolution used is 50 x 50. For the Hall 
drift modes, the propagation velocity scales linearly with 
Bq and the gradient of 'n~^, but it is independent of the 
wavenumber of the perturbation. All these properties are 
correctly reproduced. 

4-. 1.3. The nonlinear regime and Burgers flows. 

With the two previous tests, we can conclude that the 
code is able to accurately propagate the two fundamental 
modes with the correct speeds. However, the Hall drift 
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which can be discretized by an upwind conservative method 



Figure 5: Horizontal section of the evolution of the initial configu- 
ration defined by Eq. ll27ll with Bq = 10^ and k^L = tt at t = 0,2, 
and 4 (in units of tq). The shock forms at t = 2. The classical saw- 
tooth shape developed during the evolution of the Burgers equation 
is evident. 



modes are a vahd solution only in the linear regime. Let us 
consider more carefully the evolution of the By component 
in a medium stratified in the z-direction. Assuming that 
Bx = Bz = 0, it can be readily derived that the governing 
equation reduces to 



dBy ^ , , dBy 

+ giz)By 



dt 



dx 



0, 



(26) 



where g{z) = ~^( 4^en )• This is nothing but the famil- 
iar Burgers equation in the x-direction, with a coefficient 
that depends on the z coordinate. Thus, the evolution of 
the By component of the MF is governed by independent 
Burgers-like equations with different propagation speeds 
at different heights. The solution of such equation in one 
dimension is well known and has been studied for decades. 
As a last test in Cartesian coordinates we consider the 
following initial configuration: 

B^ = 0; 

By = Bocos{kxx); 



Bz = 0: 



(27) 



on the same stratified background as in the previous sub- 
section, Eq. (|24)) . which gives simply g{z) = —/SL^/tqBq 
and allows the direct comparison to the solution of the 
Burgers equation in one dimension. 

Making use of the well-known numerical techniques ap- 
plied to the Burgers equation, in our code we give an spe- 
cial treatment to the quadratic term in By in the evolution 
equation for that component (we will proceed analogously 
for the (/3-component in spherical coordinates) . The key is- 
sue is to write the Burgers-like term in conservation form: 



dBy 

dt 



d{Bl/2) 
dx 



0, 



(28) 



dt 



Ax 



(29) 



where the problem is reduced to calculate the upwind nu- 
merical fluxes at the interfaces F^^^l"^. In this case the 
wave velocity determining the upwind direction is given by 
A = g{z)By and the flux F = g{z)By/2. Using methods 
in fltix-conservative form is particularly important when 
solving prob lems with shocks or other discontinuities (e.g. 
Tord liooi)), as a non-conservative method may give nu- 
merical solutions that look reasonable but are entirely wrong 
(incorrect propagation speed of discontinuous solutions). 
The particular treatment of the time advance of this term 
(different from the general method described in Section 
3.2) has no impact on the divergence-preserving character 
of the method, since this is the component in the symme- 
try direction of the problem (the same argument applies 
for the toroidal field in axial symmetry). 

In Fig. ([5]) we show snapshots of the evolution of the 
initial conditions ((27| with k^L ^ n, Bq ^ 10'^ , and fiL = 
0.2. It follows the typical Burgers evolution. The wave 
breaking and the formation of a shock a.t t = 2tq is clearly 
captured, as would happen in the solution of the inviscid 
Burgers equation. We must stress again that this test is 
done with zero physical resistivity, i.e., in the limit Tim 
oo, which is not reachable by spectral methods or standard 
centered difference schemes in non-conservative form. 

We should also stress that the evolution of the energy 
spectrum of a Burgers-like pr oblem leads to the scalin g 
as shown for example in iTran and Dritschell ()2010l ). 

This is exactly the expected scaling, fi rst discussed in the 

context of Hall MHD in neutron stars in lGoldreich and Reisenegger 
(jl992l) , who pointed out the analogies and differences with 
turbulent velocity fields. The scaling for the power 
spectrum has been discussed in the past in terms of a Hall 
cascade, transferring energy from larger scale to smaller 
scale modes, and there is some open debate about its in- 
terpretation as a global turbulent cascade or a local steep- 
ening of some magnetic field components. 

4-.2. 2D evolution in spherical coordinates: force-free solu- 
tions 

In spherical coordinates, one of the few analytical so- 
lutions that can be used to confront numerical results is 
the evolution of pure Ohmic dissipation modes (i.e., in the 
limit TZrn — >■ 0). If the MF satisfies V x J3 = ^_B, with con- 
stant /I, and the magnetic diffusivity rj is constant, all com- 
ponents of the MF decay exponentially as cx exp(— t/r^). 
The family of solutions satisfying the above relation is de- 
scribed by radial parts involving the Bessel spherical func- 
tions and their derivatives. We test the only solution with 
a regular behavior at the center, the dipolar I — 1 function 
of the first kind. For this solution, the MF components in 
the interior of a spherical domain of radius R are: 
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Figure 6: Evolution of the purely Ohmic modes (Eg. I30II . at t = 0.2,1,2, and 3 diffusion times (t^) for the model with /iR = 11.5 and a 
resolution of 88 X 60 equally spaced grid points. In the figure we compare the analytical (lines) and numerical (crosses) radial profiles of Br 
(left) and (right) at 6» = tt/3. 



Bo = 



— 



svnx 



Bo 

Bq [ smx 

BpXi, 
2x 



cos X COS 



COS X — X sm X sm 



X 

sin a; 



cos X sm 



X 



(30) 



where x = /ir, Xi, = fiR, and Bq is a normalization factor. 
In the hmit /i — ?> we recover the solution corresponding 
to a homogeneous vertical field Bz — > Bo/ 3. 

Given this initial condition, we follow the evolution of 
the modes on several Ohmic timescales, until the field has 
been almost completely dissipated. As boundary condi- 
tions we impose the analytical solutions for Bg and B^ 
at the central cell and the external ghost cell. Note also 
that the toroidal and poloidal components of the field are 
completely decoupled in the pure res istive limit, so that 
their evolution is totally independent (|Pons and Geppertl . 
20071) . 

In Fig. |5] we compare the evolution of the numerical 
(crosses) and analytical (solid lines) solutions of Br and 
B^ at different times, for a model with /ii? = 11.5. They 
are indistinguishable in the graphic. To quantify the de- 
viation, we have evaluated the i^-norm of the deviation 
from the analytic solution AB = [B — Ban) (normalized 
to the initial MF), defined as 



(i j) I 



(31) 



where the sum is performed over all grid cells (j,j). The 
maximum of at any time for the same angular resolu- 
tion (60 cells) and three different radial resolutions of 88, 



173, and 346 grid points is 3.7 x 
1.0 X 10~^ respectively. 



10- 



1.5 X 10-4, and 



5. Evolution of a purely toroidal field in a neutron 
star crust. 

We now consider the case of a purely toroidal MF con- 
fined to the crust Tc < r < R oi a. realistic neutron star, 
with Tr — 10.8 km and R = 11.5 km. We refer the reader 



to section 4 in iPons and Geppertl (|2007l ) and section 4 in 



Aguilera et al.l (|2008l ) for details about the neutron star 
model and the electron density profiles, which depend only 
on the radial coordinate (ne(r)). For simplicity, we fix a 
constant temperature of 10* K, which gives a magnetic 
diffusivity in the range rj e [10^^ — 10^] km^/Myr. The 
boundary conditions we impose in this case are vanish- 
ing MF at both boundaries. According to the Hall induc- 
tion equation, any initial toroidal configuration will remain 
purely toroidal during the evolution, but the shape and 
location of currents can be substantially modified. The 
initial MF is given by: 



-B. 



{R 



'(r — Tc)^ sin cos f 



(32) 



where Bo is the normalization adjusted to fix the maxi- 
mum value of the toroidal field (B^^^) in the numerical 
domain to an specific value. This corresponds to an initial 
maximum value of magnetic Reynolds number TZ™^^ » 10 
for iJ^a-'^ = 10^'' G. We use this model to highlight a few 
important issues concerning energy conservation, numeri- 
cal viscosity, and shock formation. 



5.1. Energy conservation. 

Following the definitions of Section 3.5, we check the in- 
stantaneous local conservation of energy, evaluating Eq. (ITS)) 
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Figure 7: Volume-integrated energy balance as a function of time 
for a run on a 27x60 grid. At each timestep loss of magnetic energy 
(solid line) is compensated by Joule dissipation rate (dashed), within 
the a small error (dot-dashed). The initial model is given by Eq. II32II 



0.020 



with 5™="= = IQl" G. 



in the whole volume numerical domain, during a time step 
At: 



E 



87rAi 



Ay = ^^j2Ay (33) 



where the sum is performed over the cells of the domain, 
and 5-2, are local averages inside each cell of volume 
/S.V . We have used that the Poynting flux through the 
boundaries is zero, since the MF vanishes at both the in- 
ternal boundary and the star surface. In Fig. [7] we show 
the two sides of Eq. ((33)) as a function of time (dashed and 
solid lines), for a model with B™^^ = IQ^^^ G and with a 
coarse grid (27 x 60). Even for this low resolution, the er- 
ror in the instantaneous energy balance (dot-dashed) is a 
couple of orders of magnitude less than the instantaneous 
magnetic energy loss (solid) and/or the Ohmic dissipation 
rate (dashed). 

5.2. Numerical viscosity. 

It is well known that Godunov type methods provide 
stability in the most extreme conditions at the price of in- 
troducing numerical viscosity. Reconstruction procedures, 
such as the MC method we used for our code, help decrease 
the numerical viscosity in smooth regions of the flow and 
increase the spatial order of the method. However, it is im- 
portant to check that any numerical viscosity introduced 
by the method is well below the physical resistivity, in or- 
der not to affect the results. 

In Fig. [S] we compare, for different resolutions, the rel- 
ative error on the conservation of the total energy: 



.{t) = 1 



£tot{t) 
Stotit = 0) 



(34) 



where Stotit) £b{t)+ /J Qtotdt + Stotdt. The top panel 
shows ee„(t) with fixed Ng = 60 and Nr = 14,27,57,114 
(number of points in the crust) while the bottom panel 
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Figure 8: Accumulated relative error Cen 

(t) (see Eq. lO) in the 
conservation of energy during the evolution of a purely toroidal MF 
with 5™='=' = 10^* G. We show een(i) for different resolutions. Top: 
varying Nr, with fixed Ng = 60; bottom: varying Ng, with fixed Nr. 



shows results for = 57 and varying angular resolution 
Ng = 15, 30, 60, 120. The error in the energy conservation 
grows nearly linear with time in the first 10^ yr. Assuming 
that this deviation is entirely due to the numerical viscos- 
ity of the method, we can estimate the viscous timescale 
of the numerical method by assuming that 



i(i) = 1 - CXp (-i/^-num) 



t 



(35) 



and extracting Tnum from the results. 

With a resolution of 57x60 grid points (solid line in 
Fig. [5]), the error on energy conservation is of 0.25 % at 
10"* yr, which implies Tnum ~ 4 Myr. Therefore we can es- 
timate the order of magnitude of the numerical resistivity 
as 



Ar^ 



(36) 



where Ar ~ ( — Tcore}/-^r Is the average radial size of 
a cell. For Nr = 57, we obtain rynum ~ 10~^ km^/Myr. 
Comparing this value with the typical physical resistivity 
of the crust at this temperature, 77 ^ 10^^ — 10^ km^/Myr, 
we can assess that, even with a coarse grid of 57x60 grid 
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Figure 9: Evolution of a quadrupolar toroidal magnetic field confined 
into the neutron star crust. We show snapshots at t = 0, 200, 500, 
and 4000 yr. The color scale indicates the toroidal field strength, in 
units of lO^'^ G. In the figure, the thickness of the crust has been 
stretched by a factor of 4 for clarity. 



points, the numerical viscosity is at least two orders of 
magnitude less than the physical resistivity. 

We have also studied the effect on the numerical vis- 
cosity of the Courant prefactor fee, defined in Eq. (fT6|) . 
We find that varying kc does not affect significantly the 
numerical dissipation. This effect is always much smaller 
than the dependence on the spatial resolution, as long as 
the evolution is stable (typically kc < 0.1). 

5. 3. Creation of current sheets. 

The evolution equation for a purely toroidal MF, ne- 
glecting resistivity, can be cast into the following form: 



dB^ d 
at or 



ld_ 

rde 



(37) 



where 



Ay- 

Xe = 



cos 



2'Kene r sin 9 
r^c d f 1 



Aire dr 



Her^ 



(38) 



This form explicitly shows the Burgers-like character of the 
equation. In particular, a realistic profile of ne(r) has a 
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Figure 10: Angular profiles of the same model as Fig. |9l just below 
the star surface, of (in units of 10^^ G) and Jr/c (in units of 
10l2 G/km), at four different times t = 0,200,500, and 4000 yr. 
The angular resolution used is Ng = 60. 



steep gradient, which may result in the the fast formation 
of a current sheet. The location where it forms will depend 
on the initial geometry. This hidden character will thus 
become evident when TZm ^ 1- 

In order to test this limit in spherical coordinates, we 
study the same model as in the previous section but with 
a larger field, taking 5^^=^ = lO^^ G (implying 7^™^'^ w 
1000). Fig. ini shows four snapshots of the evolution of 
such configuration, where the dominant role of the Hall 
drift term leads to the fast formation of a current sheet 
in the equatorial plane. In Fig. [TU] we show a sequence 
of angular profiles of (top) and Jr/c (bottom) close 
to the star surface. Our numerical code is able to follow 
the "shock" formation and captures the thin current sheet 
with only two grid points. 

In Fig. [TT] we show the contribution of the different 
terms in the energy balance equation as a function of 
time. In the bottom panel, the sum of the magnetic energy 
(solid) and the energy dissipated by the Joule effect (short 
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Figure 11: Instantaneous (top) and time-integrated (bottom) energy 
balance for a purely toroidal MF, with 114 X 200 crustal points. 
The variation of the magnetic energy (solid line), Joule dissipation 
rate (dashes) and dissipation rate at the shock (triple-dot-dash) are 
shown together with the total energy balance (long dashes). The 
bottom panel also shows the sum (dash-dot) of the magnetic energy 
and energy dissipated by Joule effect. 



dashes) is shown by a dash-dot hne. The sum of these two 
terms is constant and equal to the initial magnetic energy 
for smooth solutions, as shown in the previous subsection. 
We can observe that now this is very accurately satisfied 
until « 200 yr, coinciding with the formation of the shock. 
After this point, the local energy balance calculated as 
above seems to be violated. About 90% of the magnetic 
energy is dissipated in 10^ yr, but the energy dissipated 
by the Joule effect is only 13% of the initial magnetic en- 
ergy, yielding to an apparent loss of total energy of ~ 77%. 
This seemingly incorrect result has a physical explanation: 
when a current sheet forms, there is energy dissipated at 
the discontinuity not accounted for by the standard Joule 
dissipation rate (right-hand side of Eq. ([55)) ). We find the 
same result with Ng = 60, 120, or 200. Our numerical 
scheme can capture the effect of an extremely thin current 
sheet (middle panel of Fig. ITO)) even if we are dealing with 
a low resolution. However, we unavoidably underestimate 



the dissipation at the discontinuity, which affects the total 
energy balance, although the implied volume is small. 

A rigorous mathematical study of the energy dissipa- 
tion rate at the shock in one-dimensional Burgers flows, 
governed by dtu + udxU — vdxxU, shows that, in the in- 
viscid limit J/ = 0, the total energy dissipation rate goes 
as 2[u]'^/3, where [u] denotes the half- jump in the variable 
u across the shock ( Tran and Dritschel . 2010l ). By anal- 
ogy with the results obtained for one-dimensional Burgers 
flows, we propose an estimate of the energy dissipated at 
the shock across a surface Sr.: 



Qs 



1 



if the discontinuity is in the radial direction, or 



DTT 



(39) 



(40) 



if it is in the angular direction. This correction is account- 
ing for the magnetic flux lost across the surface of the 
shock. The coefhcients A^, Xe are defined in Eqs. ([551) . 
Note that the correction has to be applied only at the in- 
terfaces where Aa[i3,^] is negative, which is the necessary 
condition for having a shock. In Fig. [TT]we show the the 
effect of the above correction in the energy balance when 
a discontinuity in the toroidal field is detected. Initially, 
when the MF is smooth, the extra term (triple-dot-dash) 
plays no role. It should vanish, but discretization errors 
result in a < 1% correction. However, when the current 
sheet is formed, it becomes the dominant contribution to 
explain the loss of magnetic energy. Applying this correc- 
tion, the total energy balance is very well satisfied. After 
10^ yr, when the MF has been almost completely dissi- 
pated, the error in the total energy conservation is only 

Cen - 7 X 10-4. 

Therefore, this estimation of the dissipated energy at 
the shock turns out to be an excellent approximation. Note 
that in this particular case (a purely toroidal field form- 
ing a current sheet exactly at the equator, with the dis- 
continuity in the angular direction) the generalization of 
the planar exact result to spherical coordinates is a very 
good approximation, since each radial layer is governed by 
an independent one-dimensional Burgers-like equation in 
the vicinity of the equator. In more general cases (asym- 
metric with respect to the equator, or including poloidal 
components) this approximation must be taken only as 
an estimate of the energy dissipation rate, within a fac- 
tor of two. This term is of great importance for realistic 
magneto-thermal evolution models, because it strongly en- 
hances the local deposition of heat. 

6. A general example 

The code can handle any MF geometry, confined to the 
crust or extended to the center, and has been tested in all 
the limiting cases. To conclude our study, we show results 
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Figure 12: Evolution of an initially purely poloidal magnetic field. 
The figure shows snapshots at t = 0, 2, 10, and 100 kyr. Solid lines 
show the poloidal field lines and the color scale indicates the toroidal 
field strength, in units of 10^^ G. In the figure, the thickness of the 
crust has been stretched by a factor of 4 for clarity. 



for a general case consisting in a crustal confined magnetic 
field with both toroidal and poloidal components. 

We apply vacuum boundary conditions by setting = 
for r > R and matching the poloidal field with a mul- 
tipolar vacuum solution. The vacuum solution is imposed 
as a von Neumann boundary condition problem for a po- 
tential magnetic field, expressed a.s B — V^*, where Vl/ is 
the magnetic scal ar potent ial. Using the Green's repre- 
sentation formula ( JacksonL Il991 ) with \1/ and the choice 
of Green's function G{f,r') = \f— r'\~^, the relation be- 
tween Br and Bg at surface r — R can be expressed as (M. 
Reinhardt, priv. comm.): 



/>6 fix /"TT 

47r / Bg{R,0')d0' + Be{R,0') f{e,e")d0"de' ^ 
Jo Jo Je' 

= -2 f Br{R,0')f{9,9')d0' , (41) 
Jo 



where 



n/2 



V8sin6l' d(t)' 



a/ l-cos(0- 6*') +2 sin 6* sin 61' sin'' 



Alternatively, the vacuum solution can be expressed in 
terms of Legendre polynomials (Pi) as follows: 



Br^^bi{l + l)Pi(cos( 



R 



-{1+2) 



dPi{cos0 ) (R 



-((+2) 



(42) 



where bi are the weights of the multipoles, obtained from 
the Legendre decomposition of Br{R, 9): 



21 + 1 



Br{R,0)Pi{cos 



(43) 



At each time step, we cast Br{R, 0) into Eq. (gT]) or Eq. 
to reconstruct the value of Bq in the first external ghost 
cell. We have verified that, for a given Br{R,9), both 
methods provide the same function Bg(R,0), reproduc- 
ing correctly the analytical cases (e.g., a dipole). In any 
case, the vacuum boundary condition is equivalent to avoid 
current escaping (entering) from (into) the star. The non- 
vanishing Poynting fiux across the outer boundary allows 
the (small) interchange of magnetic energy with the exter- 
nal field. For the inner boundary, we set superconducting 
boundary conditions (i.e. Br = and vanishing tangential 
electric field Eg — = 0) at the crust/core interface. As 
a consequence, the Poynting flux is zero and no energy is 
allowed to flow into/from the superconductive core. 

Our initial model consists of a purely dipolar poloidal 
component built from the following potential vector 



A^{r,9) = -f{r) sint 



(44) 



where f{r) is a radial function smoothly matching with the 
vacuum solution at r = R and vanishing r adial component 
atr = ra (sec §2.1 of lAguilera et all (|2008l) for details). We 
normalize the MF intensity at the pole to Bp — 10^^ G and 
fix the temperature to T = 10^ K during the evolution. 
This implies initial maximum values of TZrn ^ 150. 

Fig. [12] shows the initial configuration (top left panel) 
and three snapshots of the evolution of the described model. 
Because of the Hall term, the initially poloidal field imme- 
diately develops a toroidal component of opposed parity 
(mainly quadrupolar) on a timescale of about 2 x 10^ yr. 
Note that in this case the MF components must remain 
(anti)symmetric with respect to the equator. When the 
toroidal field grows to values similar to the poloidal field 
strength, the Hall term also leads to the displacement of 
MF lines to the equator and the formation of a discon- 
tinuity in 0-direction, as discussed in previous sections, 
but now in both Br and B^. After that, whistler waves 
are launched towards the poles and some sort of oscillatory 
mode appears, with transfer of energy between the poloidal 
and toroidal components. The structures are well resolved 
even with a coarse grid of 60 x 27 grid, and the disconti- 
nuities survive until the field is dissipated and the relative 
importance of the Hall term decreases, on timescales of 
~ 10^ yr. 



12 




Figure 13; Magnetic energy stored in the toroidal (solid line) and 
poloidal field (dashed) during the evolution. 



Fig.lT^shows the magnetic energy stored in the toroidal 
and poloidal component, normalized to the initial mag- 
netic energy. During the first 5 x 10'^ yr, the toroidal field 
gains energy from the poloidal field, until it reaches a 10% 
of the initial MF energy. Then, a quasi-periodic transfer 
of energy between poloidal and toroidal field is observed. 
The toroidal field strength exceeds locally the values of 
the poloidal field. The relatively low energy of the toroidal 
component is due to the small volume of the regions where 
its intensity is large. Both components acquire a marked 
multipolar behavior, with discontinuities located not only 
in the equator. The particular details about how much 
energy can be interchanged between poloidal and toroidal 
components strongly depend on both the geometry of the 
initial configuration and the MF strength. For weak fields 
(low TZm) when the resistive term dominates, the two com- 
ponents are nearly decoupled and the transfer of energy 
between them is negligible. In situations dominated by 
the Hall terms, it can be significant. We defer a deeper 
discussion about the astrophysical implications for future 
work. 

7. Conclusions 



astrophysical applications of interest. In this paper, we 
have focused on the detailed implementation of the meth- 
ods into a functioning computer code with emphasis in the 
importance of using the staggered grid and the conserva- 
tive formulation, and the need to consider the Burgers- 
like character of the equations. An extensive series of test 
problems in 2D, and for both Cartesian and spherical co- 
ordinates have been presented. These tests, that include 
the limiting cases of purely resistive evolution and purely 
Hall evolution (vanishing electrical resistivity), should be 
useful to others developing and testing codes for similar 
astrophysical scenarios. We have shown the importance 
of dissipation in the location of current sheets, which is 
us ually overlooked, or i ncorre ctly described. The analysis 



Tran and Dritschell (|201(3t ) shows that, using spectral 



methods, the number of modes needed to solve the prob- 
lem scales as TZm, which would require hundreds or even 
thousands of modes to tackle some typical problems with 
"magnetar conditions" . However, we have shown in sec- 
tion 4.1.3 that using upwind difference schemes to solve 
the equations in conservative form can satisfactorily deal 
with the purely inviscid limit with a coarse grid. We have 
also proposed a simple generalization of the exact solution 
found for the ID Burgers equation to estimate the energy 
dissipation rate at the location of a discontinuity, wherever 
it forms. This may have implications in the evolution of 
neutron stars that will be discussed in future works. 

We have shown that the code is stable and can follow 
to late times the evolution of the internal magnetic field in 
neutron stars with very low temperature, which had never 
been possible before. In future works, we will move beyond 
the developmental phase and we plan to use the code for a 
variety of applications, including the evolution of magnetic 
field in the crust of magnetars, the study of field submer- 
gence by accretion of matter, and magnetic flux expulsion 
form the core. We expect that our code will become a 
fundamental tool for astrophysical applications concern- 
ing neutron star magneto-thermal evolution. We are also 
confident that the detailed description of the algorithms 
provided in this paper will be useful to other researchers 
in the field for solving many problems in neutron stars. 



We have described a new code for the EMHD evolution 
of magnetic field of Neutron Stars. The code implements 
algorithms based on higher order Godunov methods with 
a discretization to evolve area averages of the face-centered 
components of the magnetic field (magnetic fluxes) . This 
combination results in a divergence- preserving scheme able 
to deal with discontinuities and follow the formation of 
current sheets. 

The primary motivation for developing this code has 
been the need to follow the evolution of the internal mag- 
netic field in neutron stars, including the Hall term, with 
an arbitrarily high magnetic Reynolds number. We also 
need to resolve magnetic field structures over a wide range 
of length scales and the formation of current sheets in the 
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